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Abstract. We investigate the long tim behavior of the following efficient second order in time 
scheme for the 2D Navier-Stokes equation in a periodic box: 

Aj ] n -I- j jU— 1 

+ V x (2V n - </> n_1 ) ■ V(2w n - w"- 1 ) - uAoj n+1 = /™+\ -AV> n = oj n . 

2k 

The scheme is a combination of a 2nd order in time backward-differentiation (BDF) and a special ex- 
plicit Adams-Bashforth treatment of the advection term. Therefore only a linear constant coefficient 
Poisson type problem needs to be solved at each time step. We prove uniform in time bounds on this 
scheme in L 2 , Hp er and H^ er provided that the time-step is sufficiently small. These time uniform 
estimates further lead to the convergence of long time statistics (stationary statistical properties) 
of the scheme to that of the NSE itself at vanishing time-step. Fully discrete schemes with either 
Galerkin Fourier or collocation Fourier spectral method are also discussed. 
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1. Introduction. It is well-known that incompressible fluid flows could be ex- 
tremely complex exhibiting seemingly random chaotic and/or turbulent behavior. 
Statistical approaches are necessary in order to describe such kind of complex be- 
havior (see for instance the treatises by Monin and Yaglom [35], Frisch [14], Foias, 
Manley, Rosa and Temam [12] , Lasota and Mackey [30], Majda and Wang [33] among 
others). If the the long time statistics of the system, i.e., the climate, is the subject 
of study, we then need to investigate the invariant measures of the system since it is 
the invariant measure (or stationary statistical solutions) that describes the long time 
statistics of the underlying dynamical system. Since most of these chaotic and/or 
turbulent systems are not amenable to analytical techniques at the present and in the 
near future, the issue of the development of numerical methods that are able to cap- 
ture the long time statistics becomes very important. Higher order efficient schemes 
are apparently preferred due to the long time integration needed. 

In this paper, we will focus on the development and analysis of an efficient second 
order two step numerical method that is able to capture the long time statistics for 
the following two dimensional Navier-Stokes system for homogeneous incompressible 
Newtonian fluids in the vorticity-streamfunction formulation (see for instance [36]) 

^+V ± V-Va;-i/Acj = /, (1.1) 

-Ai/} = w, (1.2) 

where uj denotes the vorticity, ip is the streamfunction, and / represents (given) 
external body forcing, v denotes the kinematic viscosity. For simplicity we will assume 
periodic boundary condition, i.e., the domain is a two dimensional torus T 2 = (0, 2-7r) x 
(0, 27r), and that all functions are mean zero over the torus. 

For analytical data, it is known that the solution is analytic in space (in fact 
Gevrey class regular due to Foias and Temam [13]), and hence Fourier spectral is the 
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obvious choice for spatial discretisation. As for time discretization, efficiency requires 
explicit treatment of the nonlinear term while stability calls for the implicit treatment 
of the diffusion term. Therefore we propose the following two step second order semi- 
implicit algorithm which treat the viscous term implicitly and the nonlinear advection 
term explicitly 

— — + V ± (2^" - V" -1 ) • V(2w" - u;™- 1 ) - vAu: n+1 = f n+ \ (1.3) 

-Aip j = uj j , j = n -(1,-4) 

Here k is the time step, and a;™" 1 , uj n , uj n+1 are the approximation of the vorticity at 
discrete time (n — l)fc, nk, (n + l)k respectively. The convergence of this scheme on 
any fixed time interval can be derived via standard methods (see for instance [34). 
There are many over-the-shelf efficient solvers of the (|1.3|) since it essentially reduces 
to a Poisson solver at each time step. This scheme falls into the category of the so- 
called implicit-explicit schemes (IMEX) p] [7] which combines second order backward- 
differentiation (BDF) and a special second order Adams-Bashforth treatment of the 
nonlinear term. However, we would like to point out that our scheme is different from 
those classical ones where the Adams-Bashforth treatment of the nonlinear term is 
in the form of linear multistep fashion of 2V ± ifj n ■ Vw" — V^V™ 1 ' w™^ 1 (see for 
instance Karniadakis, Israeli and Orszag [33], or Ascher, Ruuth and Wetton [I], or 
Varah [H]). The classical one is also known as extrapolated Gear's scheme. The new 
alternative treatment of the nonlinear advection term proves to be crucial in our long 
time stability analaysis. 

There is a long list of work on time and spatial discretization of the NSE and 
related dissipative systems that preserves the dissipativity in various forms (see for 
instance (TUI HH [HI 1221 H31 HH 123 123 EH E3 133 SSI ES] among many others). In 
particular, a general framework for the convergence of global attractor (in upper semi- 
continous fashion) for one step scheme was derived by Hale, Lin and Raugel [2l] , The 
issue of the design of numerical schemes that can capture long time statistical behav- 
ior was clearly identified by Sigurgeirsson and Stuart [33] ■ Related issue in the case 
of Hamiltonian system was discussed in Tupper [48] . It is discovered recently by the 
author that if the dissipativity of dissipative system is preserved appropriately, then 
one step numerical scheme would be able to capture the long time statistical property 
of the underlying dissipative system asymptotically in the sense that invariant mea- 
sures of the scheme would converge to those of the continuous in time system [53] . 
This general framework for the convergence of long time statistics has been applied 
to the infinite Prandtl number model for convection by Cheng and Wang [4] , the 2D 
Rayleigh-Benard convection by Tone and Wang [45], the 2D incompressible NSE by 
Gottlieb, Tone, Wang, Wang and Wirosoetisno [17] ■ The idea of preserving certain 
properties of PDE in numerical discretisation is a well-known theme (see for instance 
[32] and the references therein for preserving the sympletic structure for Hamiltonian 
systems, 41 for preserving dispersive relation in dispersive equations, as well as works 
cited above on preserving dissipativity for dissipative systems). 

The main purpose of this manuscript is to show, as long as the time-step is 
sufficiently small, that the second order two step (three level) scheme (|1.3[) is long 
time stable in the sense that we are able to derive uniform in time estimates in 
various Sobolev spaces for the solutions to the scheme. We emphasize that the time- 
step restriction we have is independent of the spatial discetization although it depends 
on the data. Hence this is different from the usual CFL condition. We also show that 
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the marginal distributions of the invariant measures of the scheme converge to those 
of the NSE (11.11) at vanishing time-step. This may be viewed as an improvement of 
our earlier result on the convergence of long time statistical properties for a first order 
classical efficient scheme for the 2D Navier-Stokes equations [17] . 

Multi-level numerical schemes are not discrete dynamical systems on the natural 
phase space. However, they can be naturally viewed as dynamical system on product 
space (see for instance the book by Stuart and Humphries [10], or Hill and Suli's work 
on approximating global attractor of sectorial evolution equation via linear multistep 
methods [22] among others). However, the limit of this dynamical system on product 
space is not the product of the original (NSE) system, and certain projections must 
be used in order to put it into a framework that is similar to the semigroup set-up 
(Hill and Suli 22] called it monoid). Although the general framework for convergence 
of stationary statistical properties presented in [53] may be modified to show the con- 
vergence of long time statistical properties using the monoid approach, here we use a 
more directly approach based on Liouville type equations as was done in the case of 
a first order scheme for the infinite Prandtl number model for convection, see Cheng 
and Wang [3]. The limit invariant measure should be concentrated on the diagonal 
heuristically since the numerical method converges. We have to consider marginal 
distributions of the invariant measures of the scheme as a dynamical system on the 
product space so that it is compartible to the phase space of the NSE. We are able to 
show that all marginal distribution/measure of invariant measure of the scheme con- 
verge to an invariant measure/stationary stationary statistical property of the NSE 
(|l.lj) (see Wang [51] for an application in terms of taking appropriate marginal dis- 
tribution and convergence of stationary statistical properties in the context of infinite 
Prandtl number limit within the Boussinesq model for convection). Therefore long 
time statistics of the scheme (in the sense of generalized time average for instance) 
converges to those of the 2D NSE as the time-step shrinks to zero. 

The rest of the paper is organized as follows. We demonstrate the long time 
stability (boundedness) of the solution to the scheme in L 2 , Hp er and H 2 er in section 
2. In section 3, we prove the convergence of the long time statistical properties. We 
briefly touch upon the issue of spatial discretization in section 4. Final remarks and 
conclusions are offered at the end. An appendix covers two technical lemmas that are 
used in this manuscript. 

2. Time uniform bounds for the semi-discrete scheme. We first recall the 
well-known periodic Sobolev spaces on fi = (0, 2ir) x (0, 27r) with average zero: 

0, <fi periodic with period 2ir in each direction j 

(2.D 



Hp™ is defined as the dual space of H™ with the duality induced by the L 2 inner 



product. We will use || ■ \\ '■= \J J n \ ■ \ 2 dx to denote the L 2 (f2) norm. 

The adoption of H™ is well-known (see for instance Constantin and Foias [6], 
or Temam |43| ) since this space is invariant under the Navier-Stokes dynamics (|1.1[) 
provided that the initial data and the forcing term belongs to the same space. 

2.1. Well-posedness of the scheme. We first demonstrate that the scheme 
(|1.3|) is well-posed in the space L 2 . 

Lemma 1. For f n+1 G L 2 , the scheme hl.S\) is well-posed on L 2 in the sense that 
there exists a unique solution in L 2 with L 2 data. Moreover, the following estimates 



hold: 



2IK+ 1 !! 2 + ^HVu/^H 2 < ^"ll + ll""" 1 !! + fc||r +1 ||H- 1 ) + ^(2||^ l |l + IK" 1 !!) 4 , Vn £2.2) 



2|| W l+1 || 2 + ^||Ac" +1 || 2 < (' ^M^ 1,1 + fc||r +1 | 



kC 2 

+ — -(2||w n || + ||w™- 1 ||) 2 (2||Vw™|| + UVu;"- 1 !!) 2 ^/! > 3. (2.3) 

7n particular, cj n ,uj n ~ 1 e L 2 implies u n+1 e Hp Pr ,uj n+2 £ H 2 er . 

Proof. It is easy to see that for uj n ,uj n ^ 1 6 L 2 , we have ip n ,ip n ~ 1 G i?p er . Hence 

||V ± (2V' n - V™" 1 ) ■ V(2w" - w"- 1 )!^-! < ^(211^11 + l^"- 1 !!) 2 (2.4) 

by the Wente type estimate (see proposition [3]) . Therefore, the scheme (ll.3[) which 
can be viewed as a Poisson type problem 



3w™ +] . „ii „„_i_i — i .„ , _ i, — ,„ „ „_r, 4w" — UJ 



n-1 



- vAco n+1 = f n+1 - V ± (2^ n - V"" 1 ) ' V(2w" - w"- 1 ) + — G Hp e l 

(2-5) 

possesses a unique solution in L 2 (in fact in Hp er ) and the solution depends continu- 
ously on the data. Therefore it is well-posed on L 2 . The bounds on the solution follow 
from taking the inner product of the equation above with and — Aw n+1 respec- 

tively, utilizing Cauchy-Schwarz as well as Wente type estimates (see Appendix). 
□ 

2.2. Time uniform bound in L 2 . Now we derive the long time stability of the 
scheme (|1.3[) in L 2 . 

There are two approaches to derive the uniform in time bounds. The first one 
utlizes a generalized G-stability property of the 2nd order BDF due to Hill and Suli 
[23] . The second one utlizies the original G-stability of the 2nd order BDF scheme 
[2"0] together with a novel two step discrete Gonwall type inequality (see lemma 2]). 

We first introduce the first approach based on a generalized G-stability of 2nd 
order BDF. 

2.2.1. Generalized G-norm and G-stability identity. We will utilize the 
G-stability of the 2nd order backward differentiation in an essential way. More specif- 
ically, we introduce the following positive definite matrix for fi > 



and we introduce a family of generalized G-norms on M as 



IGO) 

and the associated norms on R 2 valued functions as 



Vlko = V t G( M )V,VgR 2 , (2.7) 



Maw = [ V T G( A1 )V dx, V(x) : Q -> R 2 . (2.8) 

JQ, 



It is easy to see that |V| G ^ t ) is monotonically increasing in fi, and is an equivalent 
norm on M. 2 for € [0, 1] in the sense that there exist < C; < 1 < C u 0, such that 
for all fx e [0, 1], 

Ci\y\ 2 G{0) := Ci\W% < a|V| 2 G(/j) < |V| 2 < C„|V| 2 G(0) := C U |V| G < C U |V| 2 GW ,W £ R 2 , (2.9) 
^l|V|| G(0) := Q\\Vf G < QWVfcM < ||V||| a < C„||V|| 2 (0) := C U \\V\\ 2 G < C„||V|| 2 {fi)> W : SI -(SJB8) 

Moreover, we have the following identity due to Hill and Suli [23] (contained in the 
proof of their Lemma 6.1) which is a generalization of the G-stability of the BDF 
method (see for instance the classical book by Hairer and Wanner [20] ) 

/3 Is M 2 1 A-tr |2 1 | V |2 \ , ((1 +/i)«2 - 2t>i +W ) 2 

(- V2 -2v 1 + -v )v 2 + -v 2 = - \Vi\ G( , - — — | Vo| 



^^--IT-^T-^--^!^-— ,v 0|G( ^ T 

(2.11) 

where Vo = [wo,wi] T , Vi = [«i,W2] T € R 2 - This identity can be verified easily. 

2.2.2. Time uniform bound in L 2 . We are now ready to show that the scheme 
(|1.3[) is uniformly bounded in L 2 , provided that the time step is sufhciently small. In 
order to do so, we take the scalar product of (|1.3j) with 2kuj n+1 in L 2 and utilize the 
generalized G-stability of the 2nd order BDF scheme (|2.11j) to obtain, with fj, = ^ , 

II V II 2 1 My ||2 ||(l + ^"+ 1 -2^» + c"- 1 H 2 2 n+ „ + 

" "" G(,yfe) ~ TT^fc" " G(l/fc) 2(1 + i/fc) + 2^fc||Vu; || -vk\\u || 

+2k fe(2V™ - ^ n_1 ,2u; n - uj"-\uj n+1 ) = 2k(f n+1 ,u n+1 ), n = 1,2,- • ■ (2.12) 
where V„ = [w n , w n+1 ] T , and the trilinear term b is defined as 

7-L 



b{ip,4>,y)= / V ± ip-V<j>ipdx. (2.13) 
Jn 

Using the Cauchy-Schwarz type inequality, and the equivalent norm on Hp er that 
is determined by the L 2 norm of the gradient, , we majorize the right-hand side of 
(gZBP by 

2fc||r+ 1 || H - 1 ||V W "+ 1 || < ^||Va;"+ 1 || 2 + ^||r+ 1 || 2 ff _ 1 . (2.14) 



Utilising the Wente type estimate (|5.3p . and the skew-symmetry of the trilinear term 
b in the last two variables (|2.13|) . we bound the nonlinear term as 

2k b(V^(2i/j n - V™" 1 ), 2w™ - uj"-\ io n+1 ) 
= 2kb(2^ n - tp"' 1 , -(1 + vk)u n+1 + 2uo n - uj n -\uj n+1 ) 
= 2k b(2^ n ~ ip"' 1 , u n+l , (1 + vk)uj n+l - 2uj n + w™" 1 ) 

< 2fcC lu ||V ± (2y/ 1 - V" _1 )||Ki||Va; n+1 ||||(l + vk)u n+1 - 2uj n +oj n - 1 \\ 

< 5fcC t0 ||V„_ 1 ||||Va;" +1 ||||(l + iyfc)w" +1 -2w" +w"- 1 || 

< ~||(1 + vk)io n+1 - 2u: n +oj n - 1 \\ 2 + 2bk 2 C 2 w ||V„_ 1 || 2 ||Vw" +1 || 2 . (2.15) 
Relations (j2.12j) - (|2.15[) imply, under the assumption that vk < 1 

l|vj G( „ fc) - TT^II v «-illG(, fe) + (| - 25C 2 fc||v„_ 1 || 2 )fc||v w " +1 || 2 < h\\r +1 {gu) 



We are now able to prove the following long time global energy stability result: 
Theorem 1 (time uniform bound in L 2 ). Let uj n+1 be the solution of the nu- 
merical scheme hl.3\) and let f £ L°°(R + ; H~ l ), with \f\oo := |/|z,~(R+ ; _f/- 1 ) ■ Then 

2 I f I 

there exists Mofc = max{||Vo||G(i/fc)jPo} where po = J°° such that if the following 
time-step restriction is satisfied 

k ~ 5QC 2 w C u M 2 k > (2 - 17) 

then 

' ||V„|| < ||V n || G(l/fe) < M ofe , Vn > 0, (2.18) 



|V„ 



\G(uk) < (1 + v k) n H V °llG(«/fe) +P0 



1 (1 + I/fc) n 



Vn>0. (2.19) 



In particular, any ball in (L 2 ) 2 of radius p > po in the G(vk) norm, denoted Bqu,^ (p), 
is invariant under the scheme. 

Proof. The proof is straightforward by induction on n. Indeed, utlising the fact 
that the G norm is equivalent to classical norm (12.91) , we deduce from (|2.16l) that 

ll V «ll G( „ fc ) - 7^H V «-iIIg(Wc) + (~ - 25C 2 C„£ : ||V n _ 1 || G( ^ ) )fc||Va/>+ 1 || 2 < 
It is clear that ff27T9]t and (gHH hold for n = 0. 

Assuming that (|2.19p and (|2.18[) hold for n — 0, • • • , m, we then have 

^-25C 2 w C u k\\V m \\ 2 G{ „ k) >0 (2.21) 

and hence (|2.19|) remains valid for n = to + 1 which further implies the validity of 
QUED for u = to + 1. 

We remark here that Mofc could be replaced by a time-step fc independent quantity 

M = max{||V || G( , ) ,po},or M = max{^^,p } (2.22) 



due to the monotonicity of the generalized G-norm in p, and the assumption that we 
consider time-step k < 1. 

This completes the proof of theorem [TJ □ 

An immediate consequence of this theorem is the following absorbing property. 
Corollary 1 (absorbing property) . If the time-step is sufficiently small so that 

0<ki i^ckm = -- k °- (2 - 23) 

then 

^||V„|| 2 < \\\ n f G{vk) < 2p 2 , Vnfc > T (||Vo|| GW , l/loo) := ~ In f^^M.) . 

(2.24) 
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Proof. The corollary is an easy consequence of the time uniform bound (|2.19p . 
equivalence of the generalized G-norms (I2.9[) . and the fact that 1 + x > exp(cc/2) if 
x G (0, 1). D 

Remark 1 (alternative derivation of time uniform bound via discrete Gronwall 
type inequality with two steps). We now sketch the second approach of deriving 
uniform in time estimates based on the original G-stability of the 2nd order BDF 
scheme together with the two level generalized discrete Gronwall inequality presented 
in Lemma^4\ We first notice that after taking the inner product of (|1.3[) with uj n+1 , 
applying the classical G-stability for 2nd order BDF scheme (corresponding to (12. lip 
with fx = 0) , together with the same kind estimates utilized in the proof of Theorem 
[7] leads to 

\\Vm\\ 2 G - \\V m -i\\ 2 G + vk\\Vio m+1 \\ 2 + - 25C. 2 C u k\\ V m _! || 2 G )fc|| Vo/" +1 1| 2 < 

Multiplying (|2.25[) at m — n by A G (0,1) and add to (I2.25P with m = n + 1, we 
obtain, after utilizing the Poincare inequality, the equivalence of the G-norm (|2.9p . 
and denoting g n — ||V n || G ,£ = vCi\k,j3 = ( 1+ ^ Po 

(1 + e)g n+1 < (1 - X)g n + Xg^ 1 + fie (2.26) 
- 25C 2 CW l )fc||Vo/>+ 2 || 2 - A(~ - 25C 2 C u kg n - y )k\\ Vw" +1 || 2 

where we have used the fact that 

^fc(||V W "+ 2 || 2 + A||V W " +1 || 2 ) > A^fc||V„ +1 || 2 > QA^||V„ +1 || G . 

Now we assume that the following time-step restriction is satisfied: 

32C 2 max{ 5 n , 3 n ~\2/3}fc < v (2.27) 

It is then easy to see that if this time-step restriction is satisfied at n, then the same 
time step restriction is satisfied for all subsequent time thanks to Lemma [^] (|5.7p . 
Therefore the assumption in Lemma [7] is valid for all n provided it is valid initially. 
Hence the uniform in time estimates follows from Lemma [7] with a long time bound 
independent of the initial data as a result of (|5.8p . The same approach works in 
deriving uniform bound in H^ er and H^ er as well when combined with the techniques 
from the next section. 

A potential advantage of this alternative approach is the robustness of the inequal- 
ity 

2.3. Time uniform bound in H^ er and H^ er . For the purpose of proving 

the uniform in time H^ er estimates on the solution to the scheme (|1.3p . we simply 
multiply the scheme by — Aw n+1 and utilize the generalized G-stability of the 2nd 
oder BDF scheme to obtain 

H VV «llG(„ fe ) - 7^ll VV »-iH^ fe ) + 2^||A W » +1 || 2 - vk\\Vu^f 

-2k b(2ip n - ijj n -\2uj n - u"" 1 , Auj n+1 ) = -2fc(/™ +1 , Aw n+1 ), n = 1, 2, • • -(2.28) 

The right-hand side of (|2.28p can be majorized by 

2fc||/"+ 1 ||||A w "+ 1 || < — ||A W " +1 || 2 + — 1|/" +1 || 2 . (2.29) 
2 v 
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Utilising the Wente type estimate (|5.3[) . we bound the nonlinear term as 

2k b(2^ n - ip n ~\2uj n - u> n - 1 ,Auj n+1 ) 

< 2kC w \\A(2i(j n - V'"" 1 )||||V(2w™ - w n - 1 )||||Aw n+1 || 

< 5fcCg|V„_ 1 ||(2||Vw™ +1 || + ||V((1 + vk)iu n+1 - 2ui n + w n - 1 )||)||Aw n+1 || 

< i||V((l + uk)u n+l - 2u; n + W "- 1 )|| 2 + 25fc 2 C 2 ||V n _ 1 || 2 ||Aw n+1 || 2 + 8fcCJV„_ 1 |||| W "+ 1 ||^||A^+ 1 | 



< -||V((1 + vk)u n+1 - 2u n + w"- 1 )!! 2 + 25fc 2 C 2 HV,^!!! 2 !! Aw" +1 || 2 



^||A W »+ 1 || 2 + ^||V„_ 1 || 4 || W »+ 1 || 2 . 



(2.30) 



Combining the inequalities above, under the assumption that v k < 1, we obtain 



1 



I W„|| 2 - ——HVV^II 2 + (- - 25C 2 fc||V^ 1 || 2 )fc||Aa/ l+1 || 2 



'G(Wc) 1 + vk \ 

<-fciir +1 ii 2 + ^iiv„_ 1 || 4 i| w «+ i ii 2 . 



(2.31) 



Similar to the derivation of theorem [TJ we are able to show the following time 
uniform estimate in Hp er and H 2 er : 

Theorem 2 (time uniform bound in Hi er and H 2 ). Let uj n+1 be the solu- 
tion of the numerical scheme hl.3\) and let f e L°°(R + ; Hp er ), with \\f\\oc '■= 

II/IIl°=(r+ ; l 2 )' ll v /l|oc := l|V/|| L oo (R+; i2). Then there exist constants /9i(|| V ||, po, v, ||/||oo), 02(||Vo||, po, v, ||V/||oo)| 
such that if the time-step restriction (|2 . 1 T|) is satisfied, we have, 



\\W n f G{vk) < - vk)n _ 2 ||VV 2 || 2 GHfc) + Pi ( 1 1 Vo 1 1 , po i v, ll/lloc) 
l|AV n || G(yfc) < (1 + t/ 1 fc) „„ 3 |lAV 3 || 2 ;( , fc) +p 2 (||Vo||,po,^||V/|| 00 ) 



1 



1 - 



(1 + vk) n ~ 2 
1 



,n 



(1 + vk) 



n — 3 



,42^3) 



Therefore, there exist integers Ni(\\Vo\\a(vk)> Po, v, ll/lloo, k), N 2 (WVo\\ G ^ k ), Po, v, || V/||oo, fc)| 
such that 

||VV n || G(l/fe) < V2(n,Vn > Ni, ||AV„|| G( „ fe) < >/2p 2 ,Vn > N 2 . (2.34) 



Proof. The bound in H^ er is essentially there already. Indeed, if the time-step 
restriction (|2.17|) is valid, we have 

1 9 Ch 

l|vv ""°"« " TTS» vv -'»o(-) i ^ll/*' +, ll 2 + ^l|v„-,ll 4 IK +1 ll s 

(by Theorem [U and the equivalence of G-norms (|2.9|l ) 
< vkp\j2 

where 

P?:=^ll/H 2 oo + ^(l|Vo|| 6 + pg). (2.35) 
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We then deduce 



||VV n || 2 G(!/fe) <^ T ^||VV 2 ||^ fe) +p?. 



The desired estimate Hp er estimates then follows. 

As for the H 2 er bound, we multiply (|1.3[) by A 2 w" +1 and utilize the following- 
estimate on the nonlinear term: 

2k b(2ip n - iP n -\2uj n - u n -\ A 2 u n+1 ) 

< 2k\\\7 ± {2^ n - ifj n - v ) • V(2w n -w"- 1 )|| if i||VAa;" +1 || 

< 2kC w \\A{2ijj n - ^*- 1 )||||A(2w B - u;"- 1 )!! || VAw n+1 || 

< 5fc^||V„_ 1 ||(2||A w "+ 1 || + \\A{{1 + vk)u n+1 -2u n +oj n - 1 )\\)\\VAuj n+1 \\ 



1, 



< J|A((l + i/fe)w n+1 -2w" +^- 1 )|| 2 + 25fc 2 C2 || VTi _ 1 ||2|| VAci; «+i||2 + 8A;C' u ,||V„_ 1 ||||^"+ 1 ||*||VA^ ri+3 

< h A((l + uk)uj n+1 - 2uj n + w"- 1 )!! 2 + 25fc 2 C 2 ||V Il _ 1 || 2 ||VAw™+ 1 || 2 



4 



vk , — . „ Ck 
v 



— ||VAw" +1 || 2 + ^||V n _i|| 6 ||w n+1 || 2 . (2.36) 













w" 



The rest of the proof are the same as those for the Hp er estimate after we combine 
with Lemma [TJ and notice that V3 G H 2 er since V2 G Hp er . 
This completes the proof of theorem [2j □ 

3. Convergence of stationary statistical properties. The purpose of this 
section is to derive the main result of this paper, i.e., the convergence of long time 
stationary statistical properties of the scheme (| 1 . 3j) to those of the Navier-Stokes 
system (jl.ljl as time-step approaches zero. 

3.1. Dynamical system formulation. We first recall the classical approach 
of writing the two step scheme (|1.3p as a dynamical system on (L 2 ) 2 : 

(3.1) 

thanks to Lemma [1] Moreover, the rewritten scheme (13.11) can be also viewed as a 
dynamical system on Bq/ v ^\{p),\I p > po as long as the time-step restriction (|2 . 1 T[) is 
satisfied thanks to Theorem[TJ Furthermore, the trajectory is long time bounded in the 
H 1 , H 2 norm independent of the initial data thanks to the time uniform estimates in 
H 1 , H 2 that we derived in theorem[2] Therefore the dynamical system Q3.ip possesses 
a unique global attractor Ak C Bqi v ]a{pq) independent of p (see for instance Temam 
[44] ) and uniformly bounded in H 1 ,H 2 (independent of k) provided that the external 
forcing term / is time- independent. 

Due to the dissipativity of the 2nd order BDF scheme, we are able to show that 
the two components of Ak are close in the sense that the distance is no more than a 
constant multiply of the time step k. 

Lemma 2 (attractor bound and consistency) . Let f G H^ er be a time independent 
function. Then dynamical system (|3 . 1 [) is dissipative on B G ^k)(p)-:^P > Po as long 
as the time-step restriction (I2.17[) is satisfied. Moreover, the attractors are uniformly 
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bounded in the sense that for all V G Ak 
1 



1 



V|| < ||V|| G(Wc) < po, (3.2) 



VV|| < ||W|| 0(l/fc) <pi(po,Po,i/,||/||), (3-3) 



' "AV|| < ||AV|| 0(l/fc) <p 2 (p ,Po,^||V/||). (3.4) 



Furthermore, we have the following consistency estimate in the sense that there exists 
a constant Cd independent of k such that 

\\vx -v 2 \\ + y/k\\V{v x - v 2 )\\ < C d k,YV = [v 1 ,v 2 ] T G A k . (3.5) 



Proof. The uniform bound follows directly from Theorems Q] and [2j and the 
equivalence of the G-norm (j2.9[) . 

As for the distance between the two coordinates, it is straightforward from the 
scheme p. 31) , as well as the time uniform estimates derived in Theorem [5] that are 
valid on the attractor Ak, and the equivalence of the G-norm and the classical norm 

(DOS, 



- « n || < ~|K - w"- 1 !! + ^fc(||A w " +1 || + ii/n + ||v ± (2^»" - V"" 1 ) • v(2w n - w"- 1 )!!) 

< ~IK - w n_1 || + ^(||A W " +1 || + H/ll + 5<7„, || V„_i || || VV„_i||) 

< |||^" - U n ~ l \\ + |fc(||AV w || G( „ fc) + ll/H + SC^CHVn-illG^HVVn-illG^)) 

< glK ~ w n-1 || + 2 -k( P2 + ll/H + 5C w C uPoP i) 

< i-ll^ 1 - w°|| + k(p 2 + Il/H + 5C w C uPoPl ) 

provided that Vj = w J ] T , j — 0, • • • , n G Ak- Therefore 

||(Sfe Il Vo) 2 -(§fc"Vo)i|| < ^Wu 1 - l>°\\ + k(p 2 + H/ll +SC w C u p 0Pl ) (3.6) 



On the other hand, thanks to the invariance of Ak under (13.1[) . for any V G Ak 
and any integer n, there always exists a Vo G Ak so that §fc™Vo = V. Letting 
n approach infinity in the last inequality above, we deduce the desired estimate on 
Vi — v 2 . 

The desired H 1 estimate follows from interpolating the time uniform H 2 estimates 
on the solution derived in Theorem[5]and the L 2 consistency estimate obtained above. 
This ends the proof of the lemma. 

□ 

We remark that the distance between the two coordinates in the H^ er norm can 
be improved to order k as well. This desired estimate in H 1 follows from a similar 
argument by multiplying the scheme (jl.3|) by —A(uj n+1 — ui n ) and perform the usual 
energy estimates. We leave the detail to the interested reader. 
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3.2. Convergence of stationary statistical properties. We first recall the 
notion of invariant measure and stationary statistical solution that characterizes the 
long time statistics of any given dynamical system. 

We recall the definition of invariant measures. 

Definition 3.1 (Invariant measures). A Borel probability measure (j,^ on B G ^ k ^(p 
is called an invariant measure for Sk if 

[ *(S fc (V))dAi fc = f *(VW k (3.7) 

for all bounded continuous test functional ^ . 

The set of all invariant measures for Sfc is denoted TMk- 

We also recall that a Borel probability measure p on L 2 is an invariant measure, 
or stationary statistical solution, for the Navier- Stokes system if 
1. 

I ||Vw|| 2 dn(u>) < co, (3.8) 
JL 2 

2. 

( < i/Au - V^ip ■ Vw + /,$'(u) > dn(uj) = (3.9) 
Jl 2 

for any cylindrical test functional $(w) = 0((w,u>i), • • • , (ui,w m )) where (f> is 
a C°° function on M. m , {vjj,j > 1} being the standard Fourier basis which 
are also the eigenf unctions of A in L? , and <,> denotes the H , Hq per 
duality. 

3. 

[ / {v\\7u\ 2 - foj}dxdp{uj) < 0. (3.10) 
Jl 2 Jn 

The set of all stationary statistical solutions for the Navier-Stokes system is denoted 
1M. 

Roughly speaking, the first condition says that the invariant measures are sup- 
ported on the smaller and finer space of H^ er , the second condition is the differential 
form of the weak formulation of the invariance of the measure under the flow, and the 
third condition is a statistical version of the energy inequality. 

We also recall the well-known fact that the Navier-Stokes system (jl.l[) generates 
a dissipative dynamical system on L 2 with an absorbing ball of radius po (see for 
instance [44, 6]). Therefore, any invariant measure or stationary statistical solutions of 
(|l.ip must be supported in B^ipo) m L 2 . Consequently the support of any invariant 
measure or stationary statistical solution to (jl.ip must be supported on B^ 2 {po) since 
all invariant measures are supported on the global attractor (see for instance [12l[52] ). 
Henceforth, we only need to take test functionals (observables) $ to be supported on 
B L 2 (Po)- 

Next, we observe that the invariant measures of the dynamical system defined by 
our numerical scheme p.lj) and our Navier-Stokes system (|l.l[l do not share the same 
phase space. Moreover, a simple extension of the NSE to the product space is not 
appropriate, Therefore we should look at convergence of marginal distributions of the 
stationary statistical solutions to our scheme (|3.1[) . Since the attractors of the scheme 
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(|3.1[) converge to the diagonal in the product space at vanishing step-size due to the 
consistency estimates (|3.5[) , we anticipate that all marginal measures of the invariant 
measures of (|3.ip would converge to some invariant measure (stationary statistical 
solution) to the NSE (11.11) . Hence we anticipate the following result which is the 
main finding of this manuscript. 

Theorem 3 (convergence of stationary statistical properties). Let f £ Hp er be 
a time-independent function. Then the discrete dynamical systems defined via the 
scheme (|3.ip is autonomous and dissipative with non-empty set of invariant measures 
IA4k- Denote Vj,j = 1,2 the projection from (L 2 ) 2 onto its j coordinate. Let 
{/ife,fc € (0, fco]} with fif. € I.Mfc,Vfc, be an arbitrary invariant measure of the nu- 
merical scheme (|3.1[) . Then each subsequence of {fJ-k} must contain a subsubsequence 
(still denoted {/J>k}) and two invariant measure fj,j of the NSE (|l.ip so that Vj/ik 
weakly converges to fij, i.e., 



VjUk^H^k^O, (3.11) 



where V* ^ k {S) = fi k (Pr\S)),VS € B{L 2 ). 

Proof. The non-emptiness of the set of the invariant measures TMk follows from 
the uniform estimates in H 1 on the global attractors of the scheme (|3.1[) proved 
in Lemma [3.51 as well as the classical Bogliubov-Krylov technique (see for instance 

mm)- 

We observe that the family of Borel measures {fXk, k € (0, fco]} is tight in the space 
of probability measures on (L 2 ) 2 (after trivial extension to the outside of the absorbing 
ball) (see for instance Billingsley [2., or Lax [31]). Therefore, for any subsequence of 
{[J.k,k € (0, fco]}, there exists a subsubsequence, still denoted {^fc,fc £ (0, fco]}, and 
a Borel probability measure jj,o on (L 2 ) 2 such that Hk weakly converges to Our 
goal here is to show that the marginal measures of /io are in fact stationary statistical 
solutions of the Navier- Stokes system (11.11) . i.e. V*^q € TM. We work on the case 
of j = 1 without loss of generality. The proof is similar to those presented in Cheng 
and Wang [4] by utilizing the differential form of the invariant measure (stationary 
statistical solution). 

The first condition in the definition is easily verified since the global attractors 
for the discrete dynamical systems Q3.ip are uniformly bounded in Hp er independent 
of the time step fc thanks to Lemma 13.51 and the fact that the invariant measures are 
supported on the global attractor [321 [52] ■ 

In order to check the second condition, i.e., the differential form of the weak 
formulation of invariance, we let $(w) = <fi((u>, w±), • • ■ , (u>, w rn )) — </>(zi, • • ■ , z m ) be 
a smooth cylindrical test functional. Notice that 




(3.12) 



hence, denoting <, > the duality between H 1 and H~ er , we have 



I L' 2 



f <i/A«i-V- L Vi-Vwi + /,* / (wi)) rf M o(V) (V = [wi,« a ] T ,-A^ J = i; i ,i = l,2) 
J(i, 2 ) 2 

(by definition of marginal measure) 

- — / (wiAwj + V^Vi • VvijVi + fiVj) dxdfj, (V) 

(by the choice of the cylindrical test functional $ and integrations by parts) 
deb f 

2_ \ q — / (wiAtOj + V^Vi ' VwjUi + /wj) ebcd/Zfc(V) 

,l 2 ) 2 "^j Jn 



(by the definition of weak convergence) 



= lim / £ IT- I (^(V) 2 Au., + V ± (2^ 2 - ' Vwj(2v 2 - «i) + fw 3 ) dxd Mfe (V) 
fe ^° J(i, 2 ) 2 J— J 0Sj Jn 

(by the consistency estimate p.5[) as well as the choice of the cylindrical test functional) 



fe^O 



lim / (iyAS k {V) 2 -\/ ± (2^ 2 -^ 1 )-V(2v 2 -v 1 ) + f,^{v 1 )) dfi k (V) 



(L 2 ) 2 



(by the choice of the cylindrical test functional) 

= lim / /i/AS fc (V) a - V ± (2^ 2 - Vi) • V(2« 2 - m) + /, $'(^2 - ^fe(V) 
fc ^°J(L 2 ) 2 \ 2 2/ 

(by the consistency estimate p. 51) as well as the choice of the cylindrical test functional) 
(according to the scheme (|1.3[) . p. II) ) 



lim / I($( Sfc (V).[-i, |n_$(v.[-i, f] T ))^ fe (V) 
fc->-0 Jtjjivi k 2 2 2 2 



(by calculus and Lemma I3.5P 


1 3 

(by the invariance of under §/. applied to the test functional <fr(V ■ [— — , ^l 7 )) 



where we have used the boundedness and continuity of J^- on the union of the support 
of /Xfc, the consistency estimate p.5p . the invariance of /z^ under 

This proves the differential form of the weak invariance of Vi/j-q under the 2D 
Navier-Stokes dynamics, i.e., no. 2. 

The energy type inequality no. 3 can be verified easily as well. For this purpose, 
we first show that any invariant measure /ifc of the numerical scheme (|1.3|) must satisfy 
an energy type estimate. The desired continuous one will be the limit as the time 
step approaches zero. 
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Taking the inner product of the scheme (|1.3[) with u n+1 and we have 



i(||Vl 2 G - HV- 1 ^ + - 2u n - w "- l || 2 ) + H|v^ +1 | 

+ / (V ± (2V'™ - V 1 ™" 1 ) • V(2u/" - w n_1 )a; ,l+1 - /a/ l+1 ) = 0. 
in 



This can be re-written using the discrete dynamical system notation §k as 



-(||s fe (v)|| 2 G - ||V|| G + ||(s fe (v) - v) • [i, -i] T || 2 ) + H|v(§ fe (v)) 2 || 2 

+ f (V x (2^> 2 - ^) • V{2v 2 - Vl ) - /)(S fe (V)) 2 = 0. 
Jn 



Integrating this identity with respect to the invariant measure /Xfe and letting 
k — > we have 



0>-liminf-L f ||(§ fe (V)-V).[l,-l] T || 2 d / i fe (V) 

= hminf / / (^|V(S fc (V)) 2 | 2 + (V ± (2^ 2 -^ 1 ). V(2«a - « x ) - /)(S fc (V)) 2 ) dxd/ijk(V) 
fe ^° J(L 2 ) 2 Jn 

> I I (^|Vvi| 2 + (V^Vi • Vwi - /Vi) dxd Mo (V) 

J(L 2 ) 2 JO 

= / [ {V\VU)\ 2 - fu>}dxdP^ Q (L>), 

J L 2 Jn 
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where we have utilized the invariance of under and the following estimates 
/ ||V(S fc (V)) 2 f ^.(V) = / hnr £ ^^ )h V W]? 

(since {wj} form an orthogonal basis in H per ) 



(Lebesque dominated convergence theorem) 
/■ ((S fc (V)) 2 ,A Wj ) 2 

(integration by parts) 
(invariance of /i^ under S&) 

^2 



(notation following (|3.1[1 ) 

(invariance of fik under again) 
/ ||Vui|| 2 dMfe(V), 



as well as 



/ ||V Vl || 2 ^o(V)= lim Y, I { ^TW ^o(V) 

= i im lim y / (v ^ v S )2 rf Mfc (v) 

<hminfV / (V ;i' V S )2 ^ fc (V) 

= liminf/ ||V«i|| 2 d/i fc (V), 

together with the fact that on the support of /x^ 

|6(2^2 - ^i,2« 2 - vi, (S*(V)) a ) - 

< |6(2^2 - ^1, 2w 2 - Ul) (S*(V)) 2 - «i)| + |&(2^2 - 2(«2 ~ «i),«i)| + 16(2(^2 - 

< 5CjVV|| 2 (||(S fc (V)) 2 - + ||«2 - Will) 

< O(k) 

thanks to Lemma 13.51 and the Wente type estimates Lemma [3J 
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This completes the proof of the energy type inequality (no. 3 in the definition) for 
the limit probability measure y\[i§. Therefore we conclude that the limit P*fJ.o must 
be an invariant measure of the 2D Navier-Stokes system. 
□ 

Remark 2 (attractor convergence). Since the most important part of the global 
attractor must lie in the support of at least one of the invariant measures of the 
system, the convergence of invariant measures that we derived in Theorem already 
implies the convergence of an important part of the global attractor. Convergence of 
the whole global attractor (in an upper semi- continuous fashion) can be discussed as 
well following the "monoid" framework introduce by Hill and Suli fUl/. One would 



work with the so-called "monoid" [1, l] T S(t)[—^, §] : (L 2 ) 2 -> (L 2 ) 2 . This is not a 
semi-group since it does not satisfy the requirement that it is the identity map at time 
zero. We would like to point out that the convergence of the global attractors bears 
no implication on the convergence of invariant measures since knowing the support of 
the measure provides little info on the measure itself. 

4. Fully discretized case. The analysis above can be carried over to the fully 
discretized case with either Galerkin Fourier spectral approximation or collocation 
Fourier spectral approximation (see the classical books by Gottlieb and Orszag [IB] , 
or Canuto, Hussaini, Quarteroni and Zang [3], or Peyret [36] for more on spectral 
methods as well as their applications in computational fluid dynamics.). 

4.1. Galerkin Fourier spectral approximation. This subsection is devoted 
to the long time stability of the following Galerkin Fourier spectral in space and 
BDF2-AB2 in time approximation of the two dimensional Navier-Stokes equations 



' ^ n N+ N +P J v(V^(2^-^- 1 ).V(2^~^- 1 ))-,A^+ 1 = P N (r +1 ). 

(4.1) 



where Lo N ,ij) N € Vn, Vn = {all trigonometric functions on fi with frequency in each direction at most iV},| 
Pn is defined as the orthogonal projection from L 2 (f2) onto Vn- 

Just like the semi-discrete scheme (|1.3j) . we can show that the scheme (|4.1|) is 
uniformly bounded in various spaces, provided that the time step is sufficiently small. 
More precisely, we have the following: 

Theorem 4 (uniform in time bounds on the Galerkin Fourier BDF2AB2 scheme) . 
Let lo n be the solution of the numerical scheme and let f G i°°(R + ; H^ er ). 

Then the same estimates as those stated in Theorems\l\\^hold under the same time- 
step restriction (|2.17l) . In particular, the Galerkin Fourier spectral scheme (|4.1[) is 
long time stable in the sense that the solution is uniformly bounded in time and in 
truncation wave number N using either the L 2 or H^ er or H 2 er norm. 

The proof is essentially a verbatim copy of that of the semi-discrete in time case 
and we leave it to the interested reader. 

4.2. Fourier collocation approximation. We now consider Fourier colloca- 
tion spectral approximation of the semi-discrete scheme (|1.3|) . The collocation spectral 
approximation may be desirable since fast Fourier transform can be utilized to evalu- 
ate the nonlinear advection term efficiently. Here we follow a recent work by Gottlieb, 
Tone, Wang, Wang and Wirosoetisno [T7] . 

4.2.1. Fourier collocation notation. In order to present our fully discrete 
scheme with collocation spectral approximation, we need to introduce various discrete 
differential operators on these Fourier collocation spectral space. 
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For a fixed positive integer N = N x = N y , we define the mesh size h = ^ = 
h x = h y , together with the numerical grid points [xi,yj), with Xi = ih, yj — jh, 
0<i,j<N. 

Let / be a periodic function over the given 2-D numerical grid with its discrete 
Fourier expansion given by 

[AT/2] 

fi,j = ^2 fk u h exp(kiixi)exp(hiyj) , (4.2) 

k u h=-[N/2] 

Then its collocation Fourier spectral approximations to first and second order partial 
derivatives are given by 

[N/2] 

(T>Nxf) it j= E exp(i(kiXi + hyj)) , (4.3) 

fci,li=-[iV/2] 
[JV/2] 

{^nJ)^ = E (-k^fk^expiii^Xi + hyj)), (4.4) 

fci,ii=-[iV/2] 

with T>Ny , ~D % y defined analogously. In turn, the discrete Laplacian, gradient and 
divergence operators are defined as 

A N f=(V 2 Nx +V 2 Ny )f, Vjy/^ 

at the point-wise level. It is also straightforward to verify that 

V N ■ Vjv/ = A N f. (4.6) 

We also recall, for any given periodic grid functions / and g over the 2-D numerical 
grid, the spectral approximations to the L 2 inner product and L 2 norm are defined 
as 

JV-l 

II/II2 = V / (7J), with (f,g) =h 2 Y, hi9i,r (4-7) 

i,j=0 

Meanwhile, such a discrete L 2 inner product can also be viewed in the Fourier space 
other than in physical space, with the help of Parseval's identity: 

[iV/2] [JV/2] 



( v Nx f \ 






\ V N yf ) 


, Vat • 





(f,9)= E fkuljkuh = E 9ki,hfki,h> ( 4 - 8 ) 
ki,h=-[N/2] ki,h=-[N/2] 

in which fk lt ix, 9ki,h are the Fourier coefficients of the grid functions / and g in the 
expansion as in (|4.2j) . Furthermore, a detailed calculation shows that the following 
formulas of integration by parts are also valid at the discrete level: 

/> Vjv • ( f ) ) = - / Vat/, ( f ) ), (/, A N g) = - (V N f, V N g) . (4.9) 



92 ) I \ \ 92 

Discrete Sobolev spaces (with fractional order) on the 2D numerical grid can be 
defined utilizing the same idea as in the continuous case with the help of the discrete 
Laplacian operator. 
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It is well-known that the existence of aliasing error in the nonlinear term poses 
a serious challenge in the numerical analysis of Fourier collocation spectral scheme. 
Next, we recall a periodic extension of a grid function and a Fourier collocation 
interpolation operator is introduced (see for instance |17)V 

Definition 4.1. For any periodic grid function f defined over a uniform 2-D 
numerical grid, we denote f N its periodic extension. More specifically, assume that 
the grid function f has a discrete Fourier expansion in the form of 

[N/2] 

fi,j = ^2 huh exp (i(fciafj + hyj)) , (4.10) 

k u h = -[N/2] 

its periodic extension (projection) into P N is then defined as 

[N/2] 

f N (x)= /iWi «q> + litf)) • (4-11) 

k lt h=-[N/2] 

Moreover, for any periodic continuous function f, which may contain larger wave 
length, we define its collocation interpolation operator as 

[N/2] 

f i,j = (/c)fel,*i ex P (K k l X i + hVj)) , 

k u h=-[N/2] 

[N/2] 

P c N f N (x)= J2 (fc) kl ,hexv(Kkix + hy)). (4.12) 

ki,h=-[N/2] 

Note that f c may not be the Fourier coefficients of f , due to the aliasing errors. 

4.2.2. The collocation Fourier spectral BDF2 AB2 scheme. We are now 

ready to present a collocation Fourier spectral approximation in space of the 2nd order 
BDF2AB2 scheme (O]) as follows: 

— 2k N N + ^ N (m - Vn 1 ) • v« - c^- 1 ) + v w • (v£« - r N - l ){^N - ^ n ~ 1 ))) 

= v^ n + l +f™ + \ (4.13) 
-A n ^' =J N ,j = n- l,n. (4.14) 

Note that the nonlinear term is an explicit 2nd order in time spectral approxima- 
tion to alternative formulation of the nonlinear advection term as \ (V ± if!Vuj + V ■ (V^t/iw))! 
at time step n+l. This alternative formulation of the nonlinear term is due to Temam 

A straightforward application of integration by parts formula (|4.9I) gives 

(uj, u-\7 n uj + ■ (uuj)) — (uj, u-\7 n uj) — (Vjv^, uui) = 0. (4-15) 

This orthogonality property is important in deriving uniform in time and in mesh 
size estimates on the solutions to the scheme ()4.13|) . Indeed, this can be combined 
with a inductive and bootstrap argument starting with the assumption that there 
exists a constant C$ such that Hw^H^-a < Cg, for some 5 > at time step j = n— l,n 
to deduce uniform in time bounds. Detail for a first order in time scheme that treat 
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the viscous term implicitly while the advection term explicitly can be found in the 
recent work by Gottlieb, Tone, Wang, Wang and Wirosoetisno [T7]- We leave out 
the detail for the verification of the 2nd order in time scheme (14. 13)) due to space 
limitation as a future work. 

5. Conclusion and remarks. We have shown that the second order in time 
scheme (|1.3|) is long time global stable in the sense that the solutions remain bounded 
in L 2 and are asymptotically bounded in H^ er and H 2 provided that the time step 
is sufficiently small. This is a very efficient scheme since only a Poisson type solver 
is needed at each time step. This scheme is advantages over the fully explicit ones 
in terms of stability, and over the fully implicit one in terms of efficiency and unique 
solvability (the unique solvability of fully implicit scheme requires small time-step as 
well). We have also demonstrated that this scheme can be viewed as a dissipative 
discrete in time dynamical system on a product space as defined in p. II) . Moreover, 
the marginal distributions of the invariant measures of the discrete dynamical system 
(13. If) converge to invariant measures of the Navier-Stokes system (jl.ip at vanishing 
time-step. Hence the long time statistics of the numerical scheme (|1.3I) converge 
to those of the original 2D NSE (jTTTJ) . Globally long time stable fully discretized 
versions via spectral Galerkin or spectral collocation are also discussed. Therefore, 
the scheme (jl.3l) that we proposed looks very promising as a tool in the numerical 
study of long time statistical properties of the 2D Navier-Stokes system in a periodic 
box. Numerical experiments are under way to test physically interesting cases such 
as enstrophy transfer etc. 

However, there are many questions that remain to be answered. 

First, it is well-known that the set of invariant measures to the NSE (|1.1[) contains 
more than one point at high enough generalized Grashoff number in general situation 
thanks to classical bifurcation analysis result (since there are multiple steady states in 
generic case, see for instance [43]). Therefore the convergence that we derived is only 
an upper semi-continuity result. On the other hand, it is generally believed that the 
physically relevant and observable equilibrium statistical property of the 2D NSE at 
large enough Reynolds number is unique and consequently they can be characterized 
via Kraichnan's law etc (see for instance P2J HH [351 [5] among others). Hence the 
question of if the numerical scheme that we proposed above (|1.3[) can capture this 
"physically relevant one" becomes an important issue. One could argue that the 
physically relevant one should be robust under small perturbation and hence any 
reasonable numerical scheme (such as the one that we proposed in this paper) should 
be able to capture the robust behavior asymptotically. Yet another approach is to 
use small random perturbation to bring out the generic behavior of the underlying 
system since it is known that many dissipative system with appropriate noise would 
possess a unique invariant measure (see for instance Da Prato and Zabczyk [8], E 
[9], Hairer and Mattingly [19] among many others). In this case we then need to 
study numerical schemes that are able to capture the long time statistical properties 
of infinite dimensional (since we work on PDE) random (for the noise) dynamical 
systems. Besides issues that are parallel to the deterministic case, new questions 
arises such as the form of the noise as well as the magnitude etc. Of course random 
dynamical systems emerge in many other contexts including the important application 
of modeling model errors. We refrain from surveying work in this area. 

Second, the global in time stability result that is rigorously proved here imposes a 
time-step restriction (12.171) although it is of the same order (in terms of the dependence 
on the viscosity) as that for the first order scheme that we investigated earlier [IT] . We 
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believe that some kind of time-step restriction is needed due to the explicit treatment 
of the advection term. It would be interesting to find out the best possible time-step 
restriction for small v (we need to set / to be of the order of v as well if we are 
interested in having a meaningful small v asymptotics) . In practice, an on-the-fly 
criterion may be used to ensure the boundedness of the solution within the discrete 
L 2 norm for instance. 

Third, we have investigated the periodic boundary and one particular efficient 
second order in time discretization only. It would be interesting to study the case 
with the physically more interesting no-slip boundary condition, as well as other effi- 
cient implicit-explicit (IMEX) schemes such as treating the advection term explicitly 
using classical second order linear multistep approach, Crank-Nicolson leap frog type 
approach, other spatial discretization, as well as higher order in time schemes. Of 
course it is also of interest to investigate applicability of such schemes to other sys- 
tems. We note that globally long time stable schemes do exist according to Gottlieb 
and Wang [18]. However, it is not clear if those schemes are able to capture long time 
statistical properties since the modified energy used to show their stability depends 
on time-step, and the norm is not equivalent to the standard one. Higher order BDF 
may not be desirable due to the lack of A-stability [20] . 

Appendix A: Wente type estimates. We recall here a few Wente type esti- 
mates from [17] that are applicable to our doubly periodic setting. Original estimate 
of the Jacobian term (essentially H~ Y norm) goes back to Wente [54]. L 2 norm of the 
Jacobian in the case with homogeneous Dirichlet boundary condition can be found in 
Kozono and Taniuchi [35], as well as Kim [27] . 

Lemma 3. There exists an absolute constant C w > 1 such that 
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Appendix B: A discrete Gronwall type inequality for two step itera- 
tions. Lemma 4. Let {<?"} be a non-negative sequence. Suppose there exist constants 
e > 0, P > 0, A G (0, 1) such that 



5 »+ 1 <-^ 5 " + i-V- 1 + -^,Vn>l. (5-6) 
1 + £ 1 + e 1 + £ 



Then we have, for 7 = 1 ] ^^ 2 < 1, and n > 2, 



g n+1 < 7 max{ 5 ", 5 "- 1 ,2/3}, (5.7) 
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n+1 < 7 max{ 7 L: V J 5 2 , 7 LV-V ) 2/3}. (5.8) 



where \_-\ denotes the floor function (the biggest integer bounded by ■). 

Proof. The proof is a straightforward induction based on if n is even or odd and 
we leave it to the interested reader. 
□ 
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